// Copyright 1998-2016 Glenn McIntosh
// licensed under the GNU General Public Licence version 3

/** @mainpage Mathematical function library.
	@version 1.0
	@author Glenn McIntosh
	@date 1998-2016

	These classes and functions are intended for numerical applications. They
	are not complete, and are added to as required. The library uses STL
	containers, and exception handling for errors. Some functions have been
	rewritten from Pascal and C equivalents, so they are not yet template or
	class based.
	*/

// include files
#include "linear.h"
#include "matrix.h"
#include "kalman.h"
#include "quaternion.h"
#include "transform.h"
#include "smooth.h"
#include "crc.h"
#include "zeta.h"
#include "random.h"
#include "lsq.h"
#include "filter.h"
#include "ntt.h"
#include "natural.h"
#include "erasure.h"
#include <iostream>
#include <iomanip>
#include <chrono>
#include <random>
#include <complex>
template<typename T> int isize(T t) {return static_cast<int>(t.size());}

// output
template<size_t N, size_t M> void display(math::Matrix<N, M> x)
{
	for (int i = 0; i < isize(x); i++)
	{
		for (int j = 0; j < isize(x[i]); j++)
			std::cout << x[i][j] << ' ';
		std::cout << std::endl;
	}
}
void display(const math::Vector &x)
{
	for (int i = 0; i < isize(x); i++)
		std::cout << x[i] << ' ';
	std::cout << std::endl;
}
void display(const math::Vector2D &x)
{
	for (int i = 0; i < isize(x); i++)
	{
		for (int j = 0; j < isize(x[i]); j++)
			std::cout << x[i][j] << ' ';
		std::cout << std::endl;
	}
}
void display(const math::Mpint &x)
{
	for (int i = 0; i < isize(x); i++)
		std::cout << x[i] << ' ';
	std::cout << std::endl;
}

// test routine
void test()
{
	try
	{
		math::Vector a(128, 0.0);

		// initial data
		math::UniformRandom x(1);
		for (int i = 0; i < isize(a); ++i)
			a[i] = float{x} * 100;
		std::cout << "*** vector data" << std::endl;
		display(a);

		// test correlation
		std::cout << "*** correlation" << std::endl;
		math::Vector c = a;
		math::autoCorrelate(c);
		display(c);

		// test wavelet transform
		std::cout << "*** Haar wavelet transform data" << std::endl;
		math::haarTransform0(a);
		display(a);

		// test inverse wavelet transform
		std::cout << "*** inverse Haar wavelet transform data" << std::endl;
		math::haarTransform1(a);
		display(a);

		// test hartley transform
		std::cout << "*** Hartley transform data" << std::endl;
		math::hartley(a);
		display(a);

		// test inverse hartley transform
		std::cout << "*** inverse Hartley transform data" << std::endl;
		math::hartley(a);
		display(a);

		// test smooth
		std::cout << "*** polynomial smoothed data" << std::endl;
		math::Vector b = a;
		auto fir = math::savitsky_golay_filter<10, 2>(0);
		math::Vector response{fir.begin(), fir.end()};
		math::fir_smooth(b, response, true);
		display(b);
		std::cout << "*** Gaussian smoothed data" << std::endl;
		b = a;
		math::Vector xx(b.size(), 0.0);
		for (int i = 0; i < isize(xx); ++i)
			xx[i] = i;
		math::gaussian_smooth(xx, b, 5.);
		display(b);

		// test least squares curve fit
		std::cout << "*** least squares curve fit" << std::endl;
		b = a;
		math::Lsq<2> lsq;
		lsq.add(xx, b);
		lsq.fit();
		lsq.evaluate(xx);
		display(xx);

		// test CRC
		std::cout << "*** CRC ***" << std::endl;
		const char *test0 = "CatMouse987654321";
		math::Crc t0;
		std::cout << std::hex;
		std::cout << "CRC=" << (t0 += (const uint8_t *) test0) << std::endl;

		// test zeta function
		std::cout << "*** Zeta ***" << std::endl;
		Zeta<24> zeta;
		std::cout << zeta(0.0) << std::endl;
		std::cout << zeta(0.5) << std::endl;
		std::cout << zeta(1.0) << std::endl;
		std::cout << zeta(1.5) << std::endl;
		std::cout << zeta(8.0) << std::endl;
		std::cout << zeta(-1.0) << std::endl;
		std::cout << zeta(-2.0) << std::endl;
		std::cout << zeta(-7.0) << std::endl;
		std::cout << zeta(-8.0) << std::endl;
		std::cout << zeta(std::complex{0.5, 0.5}) << std::endl;
		std::cout << zeta(std::complex{0.5, 13.0}) << std::endl;
		std::cout << zeta(std::complex{0.5, -21.022040}) << std::endl;
		std::cout << zeta(std::complex{0.5, -25.010858}) << std::endl;

		// test random number generators
		std::cout << "*** hashing and random number generators ***" << std::endl;
		math::UniformRandom r0(1);
		float{r0};
		unsigned long int h = r0.Hash(1);
		std::cout << h << ' ' << float{r0} << std::endl;
		h = r0.Hash(99);
		r0.SetIndex(99);
		std::cout << h << ' ' << float{r0} << std::endl;
		math::UniformRandom r1(99);
		float{r1};
		h = r1.Hash(1);
		std::cout << h << ' ' << float{r1} << std::endl;
		h = r1.Hash(99);
		r1.SetIndex(99);
		std::cout << h << ' ' << float{r1} << std::endl;
		math::RandomBitSequence r2(0xBAA96887);
		for (int i = 0; i < 64; ++i)
			std::cout << (r2?'1':'0');
		std::cout << std::endl;
		std::cout << std::dec;

		// test filter
		std::cout << "*** IIR filtering ***" << std::endl;
		math::IirFilter f0(4);
		f0 = 345;
		for (int i = 0; i < 30; ++i)
			std::cout << (f0 += 346) << ' ';
		std::cout << std::endl;

		// initial data
		math::Matrix<3,3> m{};
		for (int i = 0; i < isize(m); ++i)
			for (int j = 0; j < isize(m[0]); ++j)
				m[i][j] = float{x} * 100;
		std::cout << "*** matrix data" << std::endl;
		display(m);

		// invert data
		std::cout << "*** multiplied by inversion" << std::endl;
		display(math::operator*(m, math::invert(m)));

		// normalize rotation matrix
		std::cout << "*** normalize rotation matrix" << std::endl;
		display(math::normalizeRotationMatrix(m));

		// square of cholesky decomposition
		{
			// initial data
			math::Matrix<4,4> m{};
			m[0][0] =  5.0; m[0][1] =  1.2; m[0][2] =  0.3; m[0][3] = -0.6;
			m[1][0] =  1.2; m[1][1] =  6.0; m[1][2] = -0.4; m[1][3] =  0.9;
			m[2][0] =  0.3; m[2][1] = -0.4; m[2][2] =  8.0; m[2][3] =  1.7;
			m[3][0] = -0.6; m[3][1] =  0.9; m[3][2] =  1.7; m[3][3] = 10.0;
			std::cout << "*** matrix data" << std::endl;
			display(m);

			// cholesky decomposition
			math::choleskyDecompose(m);
			std::cout << "*** cholesky decomposition" << std::endl;
			display(math::operator*(m, math::transpose(m)));
		}

		// test number theoretic transform
		{
			math::Mpint a(128);

			// initial data
			math::UniformRandom x(1);
			for (int i = 0; i < isize(a); ++i)
				a[i] = (uint32_t) x & 256-1;
			std::cout << "*** vector data" << std::endl;
			display(a);

			// test number theoretic transform transform
			std::cout << "*** number theoretic transform" << std::endl;
			math::ntt<>(a);
			display(a);

			// test inverse number theoretic transform transform
			std::cout << "*** inverse number theoretic transform" << std::endl;
			math::intt<>(a);
			display(a);
		}

		// test natural number transformations
		{
			std::cout << "*** arbitrary precision multiplication ***" << std::endl;
			math::N x{{0xFFFF0000, 0x12345678}}, y{{0xFFFF1234, 0x12345678}};
			math::N z = x*y;
			for (int i = 0; i < isize(z.value); ++i)
				std::cout << std::hex << z.value[i] << std::dec;
			std::cout << std::endl;
		}

		// test erasure codes
		{
			// initialize data set
			constexpr int n_bytes{12};
			struct packet {uint8_t bytes[n_bytes];};
			constexpr int n_packets{16};
			uint8_t data[n_packets*n_bytes], recovered_data[n_packets*n_bytes];
			std::default_random_engine generator(std::random_device {}());
			generator.seed(1);
			std::uniform_int_distribution<> distribution(0, 255);
			for (int i = 0; i < n_packets*n_bytes; ++i)
				data[i] = distribution(generator);

			// display data set
			std::cout << "*** packet data" << std::endl;
			std::cout << std::setfill('0');
			for (int i = 0; i < n_packets; ++i)
			{
				for (int j = 0; j < n_bytes; ++j)
					std::cout << std::hex << std::setw(2) << int{data[i*n_bytes+j]} << std::dec;
				std::cout << std::endl;
			}
			std::cout << std::setfill(' ');

			// run test
			math::RemainderErasureTransmit<packet, 32> construct;
			math::RemainderErasureReceive<packet, 32> reconstruct;
			//math::LubyTransformTransmit<packet, 32> construct;
			//math::LubyTransformReceive<packet, 32> reconstruct;
			int transmit = 0;
			while (transmit < 100)
			{
				// construct packet
				packet packet = {0};
				int seed = generator()%2048;
				construct(seed, sizeof(data)/sizeof(packet), packet, &data);

				// send selection
				++transmit;
				if (reconstruct(seed, sizeof(recovered_data)/sizeof(packet), packet, &recovered_data))
					break;
			}

			// display data set
			std::cout << "*** remainder erasure code reconstructed data" << std::endl;
			std::cout << "packets=" << transmit << std::endl;
			std::cout << std::setfill('0');
			for (int i = 0; i < n_packets; ++i)
			{
				for (int j = 0; j < n_bytes; ++j)
					std::cout << std::hex << std::setw(2) << int{recovered_data[i*n_bytes+j]} << std::dec;
				std::cout << std::endl;
			}
			std::cout << std::setfill(' ');
		}

		// speed trials
		#ifdef SPEEDTRIAL

		// hartley transforms
		std::cerr << "*** hartley transforms" << std::endl;
		std::cerr << "length, MFLOPS" << std::endl;
		for (int n = 16; n <= 1*1024*1024; n *= 2)
		{
			math::Vector a(n);
			math::UniformRandom x(1);
			for (int i = 0; i < n; ++i)
				a[i] = float{x} * 100;
			int repeats = 1024*1024/n;
			if (repeats == 0) repeats = 1;
			auto t0 = std::chrono::high_resolution_clock::now();
			for (int i = 0; i < repeats; ++i)
			{
				math::hartley(a);
				math::hartley(a);
			}
			auto t1 = std::chrono::high_resolution_clock::now();
			std::cerr << n << ", " << 2.5*n*log2(n)/1e6/(std::chrono::duration<double>(t1-t0).count()/2/repeats) << std::endl;
		}

		// multiprecision multiplication
		std::cerr << "*** arbitrary precision multiplication" << std::endl;
		std::cerr << "length, MFLOPS" << std::endl;
		for (int n = 16; n <= 1*1024*1024; n *= 2)
		{
			math::Mpint a(n);
			math::UniformRandom x(1);
			for (int i = 0; i < n; ++i)
				a[i] = (uint32_t) x & 16-1;
			int repeats = 1024*1024/n;
			if (repeats == 0) repeats = 1;
			auto t0 = std::chrono::high_resolution_clock::now();
			for (int i = 0; i < repeats; ++i)
			{
				math::ntt<>(a);
				math::intt<>(a);
			}
			auto t1 = std::chrono::high_resolution_clock::now();
			std::cerr << n << ", " << 2.5*n*log2(n)/1e6/(std::chrono::duration<double>(t1-t0).count()/2/repeats) << std::endl;
		}

		// erasure codes
		{
			std::cerr << "*** Chinese remainder erasure code" << std::endl;
			// initialize data set
			constexpr int n_bytes{12};
			struct packet {uint8_t bytes[n_bytes];};
			constexpr int n_packets{16};
			uint8_t data[n_packets*n_bytes], recovered_data[n_packets*n_bytes];
			std::default_random_engine generator(std::random_device {}());
			std::uniform_int_distribution<> distribution(0, 255);
			for (int i = 0; i < n_packets*n_bytes; ++i)
				data[i] = distribution(generator);

			// test measurements
			constexpr int n_repeats{100000};
			int sum{0};
			int max{0};

			// run test
			auto t0 = std::chrono::high_resolution_clock::now();
			for (int repeat = 0; repeat < n_repeats; ++repeat)
			{
				math::RemainderErasureTransmit<packet, 32> construct;
				math::RemainderErasureReceive<packet, 32> reconstruct;
				//math::LubyTransformTransmit<packet, 32> construct;
				//math::LubyTransformReceive<packet, 32> reconstruct;
				int transmit = 0;
				while (transmit < 200)
				{
					// construct message
					packet message = {0};
					int seed = generator()%2048;
					construct(seed, sizeof(data)/sizeof(message), message, &data);

					// send selection
					++transmit;
					if (reconstruct(seed, sizeof(recovered_data)/sizeof(message), message, &recovered_data))
						break;
				}

				sum += transmit;
				if (transmit > max) max = transmit;

				// check result
				for (int i = 0; i < n_packets*n_bytes; ++i)
					if (recovered_data[i] != data[i])
					{
						std::cerr << "error: at repeat " << repeat << ' ' << transmit << std::endl;
						repeat = n_repeats;
						break;
					}
			}
			auto t1 = std::chrono::high_resolution_clock::now();
			std::cerr << "time/packet=" << 1e06 * std::chrono::duration<double>(t1-t0).count()/n_repeats/n_packets << "us" << ' ' << "average_sent=" << 1.*sum/n_repeats << ' ' << "maximum_sent=" << max << std::endl;
		}
		#endif // SPEEDTRIAL
	}
	catch (std::exception &error)
	{
		std::cerr << "exception " << error.what() << std::endl;
	}
}
// main
int main(int argc, char **argv)
{
	// flags
	int iArg = 1;
	char *arg = argv[1];
	while (iArg < argc && (arg[0] == '-' || arg[0] == '/'))
	{
		char c = tolower(arg[1]);
		switch (c)
		{
		case '?': case 'h':
		default:
			std::cerr << "Usage(V1.0): " << argv[0] << " [-options] filename" << std::endl;
			std::cerr << "   -h	 : help" << std::endl;
			return 1;
		}
		arg = argv[++iArg];
	}

	// files
	while (iArg < argc)
	{
		arg = argv[++iArg];
	}
	if (iArg < argc)
	{
		std::cerr <<"Error: too many arguments" << std::endl;
		return 1;
	}

	// test function
	test();
}
